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Abstract 

Implicit particle-in-cell codes offer advantages over their explicit counterparts 
O \ in that they suffer weaker stability constraints on the need to resolve the higher 

^ \ frequency modes of the system. This feature may prove particularly valuable for 

O ' modeling the interaction of high-intensity laser pulses with overcritical plasmas, 

c/3 ■ in the case where the electrostatic modes in the denser regions are of negligi- 

ble influence on the physical processes under study. To this goal, we have de- 
veloped the new two-dimensional electromagnetic code ELIXIRS (standing for 
ELectromagnetic Implicit X-dimensional Iterative Relativistic Solver) based on 
^ the relativistic extension of the so-called Direct Implicit Method [D. Hewett and 

^ ■ A. B. Langdon, J. Comp. Phys. 72, 121(1987)]. Dissipation-free propagation of 

light waves into vacuum is achieved by an adjustable-damping electromagnetic 
solver. In the high-density case where the Debye length is not resolved, satisfac- 
tory energy conservation is ensured by the use of high-order weight factors. In 



O \ this paper, we first present an original derivation of the electromagnetic direct 

^ ■ implicit method within a Newton iterative scheme. Its linear properties are then 

investigated through numerically solving the relation dispersions obtained for 
both light and plasma waves, accounting for finite space and time steps. Finally, 
^ \ our code is successfully benchmarked against explicit particle-in-cell simulations 

\ for two kinds of physical problems: plasma expansion into vacuum and relativis- 

tic laser-plasma interaction. In both cases, we will demonstrate the robustness 
of the implicit solver for crude discretizations, as well as the gains in efficiency 
which can be realized over standard explicit simulations. 
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1. Introduction 



Particle-in-cell (PIC) codes have become widely used plasma simulation tools 
owing to their ability to mimic real plasma behavior. Yet the standard PIC algo- 
rithm employs an explicit time-differencing, and hence suffers from strict stability 
constraints on the time step, which needs to resolve the highest-frequency modes 
of the system Furthermore, the mesh size must be comparable to the Debye 
length Ad in order to prevent the finite-grid instability [H. As a consequence, ex- 
plicit PIC codes may find it difficult to cope with the large spatial and temporal 
scales associated with a number of physical scenarios, thus requiring massively 
parallel computing facilities Several alternatives have been developed over the 
past decades to relax these constraints so that the choice of the space and time 
steps can be dictated by physical accuracy rather than stability conditions. The 
simplest way to do so is to suppress high-frequency processes within the mathe- 
matical model itself. Codes based on the Darwin-field approximation 0,01, gy- 
rokinetic equations [sl or hybrid particle-fiuid models 0, Q, S, 0, rely precisely 
on such an approach. The shortcoming inherent in these codes is the somewhat 
uncertain domain of validity of their basic assumptions. A second, more involved 
numerically, possibility retains a fully kinetic and electromagnetic description by 
using an implicit scheme for the entire Vlasov-Maxwell set of equations. This is 
the approach dealt with in this work. 

The main feature, and difficulty, of a fully implicit PIC scheme is the pre- 
diction of the future particles' charge and current densities as functions of the 
future electromagnetic fields. Two main techniques have been designed to this 
goal. The first one to be published, the so-called moment method, makes use of 
the fiuid equations to predict future source terms 111. Il2l. Il3l . 114 lisl . Ild |. and 
has been recently extended to the relativistic regime [17|. The present article 
will focus on the alternate approach, referred to as the direct iinplicit method. 



which is based on a direct linearization of the Lorentz equations [18|, ll9|, l20|, 121 



Most implementations of the direct implicit method start with the so-called D 



discretization of the Lorentz equation, first presented in Ref. [22|. The rela 



tivistic formulation, originally derived in Re f. |23l|. was implemented, albeit in a 
simplified form, in the LSP code 



24, 25, 26 




The direct implicit method proceeds as follows. First, particles' momenta and 
positions are advanced to an intermediate time level using known fields, yielding 
predicted charge and current densities. Second, by linearizing the latter quantities 
around the predicted momenta and positions, we can express correction terms as 
functions of the future fields and thus derive an implicit wave equation. Once 
this equation is solved, the particles' quantities are updated. Here we will show 
that the direct method can be derived as a simplified Newton scheme. 

Our main motivation is the simulation of the interaction of an ultra-intense 



2 



laser pulse with solid-density plasma slabs. The energetic particle beams originat- 
ing from this interaction s tir grea t interest in many fields spanning inertial con- 
finement fusion 
nuclear physics 



mm 



38, 



39 



30, IMI, l32, l33|| , high e nerg y density physics [3J, l35|, l36|, l37 



lerg - 



or medical physics [40||. For the high plasma densities 
considered, the electron plasma frequency Up largely exceeds the laser frequency. 
Using an explicit PIC code, the space and time steps should resolve the high- 
frequency electron plasma modes of the plasma bulk. However, these modes are 
of no interest for the problem since they do not affect the laser-plasma inter- 
action nor other potentially important related processes as the subsequent, fast 
electron-driven ion expansion. By contrast, resorting to an implicit scheme would 
allow a significantly increased time step, that is, determined only by the need to 
resolve the incoming laser wave. In this respect, one should realize that the strong 
wave damping inherent with implicit methods may be harmful in the context of 
laser-plasma interaction, for which light waves have to travel over many wave- 
lengths. This prompted us to develop an electromagnetic solver with adjustable 
damping, based on a generalization of the scheme initially proposed by Friedman 



41| for the Lorentz equation. We will demonstrate that our adjustable damping 



scheme tolerates abrupt spatial jumps in the controlling parameter. Our code 
therefore allows for dissipation-free laser propagation into vacuum, along with 
strong damping of undesirable plasma waves into the densest part of the target. 

As explicit codes, implicit codes suffer from the artificial heating arising from 
a crude discretization of the Debye length, as is commonplace when handling 
large-scale, high-density plasmas. This detrimental effect is generally attributed 
to the so-called grid-instability [H. To keep it at an acceptable level, we will 



42, 43 



exploit the well-known mitigating influence of high-order weight factors 
by using quadratic weight factors. We will also take advantage of the stabilizing 
effect of the large time steps allowed by the implicit scheme. 

The paper is organized as follows. In Sec. [2], we recall the basic princi- 
ples of the PIC technique, give the implicit time-discretized equations to solve, 
and derive within a simplified Newton formalism the relativistic direct implicit 
method. In Sec. [3l we outline the numerical resolution of the wave equation 
as implemented in our newly developed, 2Dx-3Dv code ELIXIRS (ELectromag- 
netic Implicit X-dimensionnal Iterative Relativistic Solver). The introduction of 
implicit injecting/outgoing boundary conditions for the electromagnetic field is 
also discussed. Sec. [4] is devoted to the linear properties of the direct implicit 
method through the resolution of the electromagnetic and electrostatic disper- 
sion relations. The effects of finite space and time steps, adjustable damping and 
high-order weight factors will be accounted for. Finally, in Sec. [5], our code is 
benchmarked against explicit simulations for two kinds of physical problems: the 
expansion of a plasma slab in vacuum, and the interaction of an ultra-intense 
laser pulse with an overcritical plasma target. The sensitivity of the simulation 
results to the damping parameter and the number of macro-particules will be 
addressed. 
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2. The relativistic direct implicit method as a simplified Newton scheme 



In contrast to Ref. 23 



we present here a derivation of the electromag- 
netic direct implicit method for the relativistic case within a Newton iterative 
scheme and a weak formulation of Maxwell's equations. Anticipating our need of 
a dissipation-free propagation of light waves inside the vacuum region of the sim- 
ulation domain, we introduce a generalization of the adjustable damping scheme 



proposed and used in the electrostatic regime by Friedman |4ll| 



2.1. Basic equations 

Consider Maxwell's equations 



V X E = - 



dB 

'dt 



V X B 



. 1 dE 



(1) 
(2) 



and the collisionless Vlasov equation for the distribution function /^(x, u, t) of 
the sth particle species 



dfs ^ udfs 
dt 7 dx. 



qs_ 



E 



u 



7 



X B 



du 



0. 



(3) 



Here and are the charge and the rest mass of the sth particle species, 
respectively, u denotes the relativistic momentum normalized by m^. The rel- 

1 /2 

ativistic factor then writes 7 = (1 + /(?) . The particle method consists in 
describing the distribution function fs as an ensemble of macro-particles in the 
form 



Ns 



/,(x, u, t) = ^ 5(x - X,(t))5(u - Up(t)) 



(4) 



p=i 



where S is the shape function Ns the total number of particles of the sth 
species, and 5 the Dirac distribution. The relativistic motion of each macro- 
particle obeys the following equations: 



dt 
d\Jp{t) 
dt 



^^E[X,(t),t] + 



Ml 



X B[Xp(t),t] 



(5) 
(6) 



We now make use of the implicit scheme with adjustable damping proposed 
by Friedman [4]J for an electrostatic problem, which generalizes the so-called Di- 
scheme of Langdon et al. [l8|, [l9|, l20|, l23| . The equations of motion are discretised 
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as 



X„ + At- 



U, 



n+l/2 



7n+l/2 

TT TT , , A N , /U„+i/2 + U„„i/2 
Lln+1/2 = U,i_i/2 + — (a„+i + A„„ij + — ( ) X tin{-^nj , 



(7) 



2m ^ 



In 



1 ^ ) , 



1 - a„ + -^SL 



f: 



a„_i . - 2 y ■ 2 
where the index n denotes the time step index and we have defined 



In 



U,i_i/2 + — (a„+i + A„_i) 



1 2 



1/2 



7n+l/2 



•-^ n+l/2 



1/2 



(8) 
(9) 

(10) 



(11) 
(12) 

(13) 



Friedman's scheme can be readily applied to Maxwell's equations, which yields 

At. 



En+l = E„ + C AtV X B„+i/2 Jn+l/2 , 

^0 

Bn+1/2 = B„_i/2 X (E„+i + E„_i) 

Bn = B„_i/2 X E„ , 

En-1 = "^E„ + f 1 ^ 1 Eri_2 



Eri-1 







E 



n-2 



(14) 
(15) 
(16) 
(17) 

(18) 



where j denotes the current density. 

As will be demonstrated in Sec. [H this scheme allows, via the parameter 
9f, a fiexible control of the damping of the high-frequency (electrostatic and 
electromagnetic) waves of the system. This property is of major interest for 
applications such as laser-plasma interaction involving a traveling electromagnetic 
wave into vacuum, for which the numerical damping associated with the standard 
Di method may prove too severe. It is worth noting that, even though referred to 
uniquely as 6f, the damping parameters involved in the electromagnetic scheme 
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and the particle pusher may assume distinct values. The next sections will be 
devoted to the solution of the set of Eqs. fl7l)- f|T8ll within a Newton iterative 
scheme. We will show that for a proper choice of the initial conditions, this 
scheme reduces to the direct implicit method developed in Refs. 



20, 23 



2.2. Weak formulation of the electric field equation 

By replacing Eq. (fTSl) into Eq. (fT4l) . one obtains the following wave equation 

c^At^ At 
E„+i + V X V X E„+i + — j„+i/2 = Q' , (19) 
z eo 

with the (known) source term 

c^At"^ 

Q' = E„ + c^AtV X B„_i/2 —V X V X E„_i . (20) 

For any test function ip, we assume the following weak formulation of the current 
density 

j„+i/2(x)?/;(x)dx 

= X] y y /s,o(x, u)V„+i/2(x, u) [tp (X„+i(x, + ^ (X„(x, u))] dxdu , (21) 
where /^ o = /s(x, u, 0) is the initial particle distribution function and 'Vn+1/2 = 

U„+l/2/7n+l/2- 

The problem then consists in finding (E„+i, X„+i, U„+i/2) which solve 

r c^At^ f 

/ E„+i(x)^(x)c/x + / V X V X E„+i(x)^(x)dx 

+ ^ J j„+i/2(x)^/'(x)dx = J Q'(x)V'(x)dx (22) 

together with Eqs. (j7l)- (fT3l) . We employ the Newton method to solve this system: 
for each quantity of interest Y, we introduce the ansatz 

yik+l) _ y{k) xyW /. _ n 1 (r,o\ 

where a = (1/2, 1) depending on whether Y is centered at full or half time steps. 
The subscript n + l will be hereafter omitted for clarity. Substituting the above 
ansatz into Eq. ( !20l) yields 

/ [E('=)(x) + 5E(^^)(x)] ^(x)c/x+ / V X V X [e('=)(x) + 5E('=)(x)] 7/;(x)rfx 

+ ^ J j^'+'Hx)V'(x)dx = J Q'(x)^(x)rfx. (24) 
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The term involving j'^'^"'"^^ is calculated with positions X'^'^+^) and velocities V^^"*"^) 
I j('=+iV(x)rfx = E I / ")^^'^ [^(^^'^) + ^(^")] ^^^^ 

+ E I / ")^^'^ [VV^lX^'^)) ■ 5XW] rfxrfu . (25) 

To obtain the equation solved for the electric field, we need to express the 
terms X'-'^-*, ^X'-'^-*, V^'^-' and (5V^'^^ as functions of the electric field. Before pro- 
ceeding, let us first define the following quantities 



R(XJ = 



1 + 



2 

' 1 + 02 
1 



1/2 



1/2 



-B„(X„) , 

i+e^e-ex i 



7 



(fc) 



I 



n(eW(x('=)),u('=)) 



qsAt 



4:171. 



Un-1/2 + U^"^ 

r(fe)3 



X B4X„) 



U„_v2 + ^f^E('=)(X('=)) + A„_i') 

/I V y 



4 

with I the identity matrix. Straightforward calculations then yield 

x<') = x. ^'u<'> 



7 



(fc) 



(26) 

(27) 

(28) 
(29) 
(30) 



(31) 



(32) 
(33) 
(34) 
(35) 



Using the above expressions and the Newton ansatz (j23ll . the Lorentz equation 
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becomes 

2 " ^ 2m, V y 

_ ^N(E('=)(X('=)), U(^))VE«(X('=))5X(^) 
2ms 

_ ^N(E('=)(X('=)), U(^))5E(^)(X(^)) , (36) 
2ms 

where we have dropped second-order terms. Assuming the electric field gradient 
term is negligible, this equation further simplifies as 



A„_i + i^E(^)(X(^)) 



m.o 



u(^=) + 5u(^) = u„_i/2 + ^ [I + R(x„,)] 

+ ^ [I + R(x„)] [I - N (E(^-)(X(^-)), U^'^))] 5E('=)(X('=)) . (37) 
4ms 

The set of equations (I22l ) - (l38l) constitutes the weak formulation of the problem. 
We will now show how to recover the direct implicit method as a simplified 
Newton algorithm. 

2.2.1. The direct implicit method 

The simplest scheme consists in considering only one iteration in the above 
system and choosing the following initial values 

= X„+i r 5X(0) = SX 

= U„+i/2 =6V (38) 

= [ 5E(0) = E(i) = E„+i , 

where we have introduced the predicted position and momentum X,j+i and V„_^i 
computed from the known fields A„_i and B„. We have 

X„+i = X„ + AtS^ , (39) 

7n+l/2 

Un+1/2 = R(X„)U„_i/2 + ^ [I + R(X„)] A„_i . (40) 
with 7n+i/2 = 7^°^- The correction terms then write 

<5U = ^[I + R(XO][I - N(U„+i/2)]E„+i(X„+i) , (41) 

6Y = M6V , (42) 
5X = AtMdV , (43) 
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where we have defined 



/i(U„+i/2) = N(0,U„+i/2 



U„_i/2 + U„+i/2 



In 



X B„(XJ 



U„_l/2 + ^^n-l ] , (44) 



and 7n = T^^K After substituting the above equations into (l25l) . using X„ 
X„_i_i — AtV„+i/2 and replacing the resulting expression into ([24l) . we obtain 



E„+i(x)V'(x)rfx + / V X V X E„+i(x)?/'(x)dx 



s 

E 



/,,o(x, u)V„+i/2(x, u) tpiXn+i (x, u)) + V^(X„ (x, u)) 



/,,o(x, u)5V(x, u)?/'(X„+i(x, u))cix(iu 



(ixdu 



Q'(x)^/^(x)dx. 



(45) 



From Eq. ( l2Ti ). we identify 



2eo 
At 



/s,o(x, U)V„+1/2(X, U) ^(X„+i (X, U)) + l/j{Xn (x, u)) 
jn+l/2(x)^/'(x)dx. 



c/xfiu 



(46) 



To reduce the next integral, it is convenient to introduce the weak formulation 
of the predicted charge density 

j ps(x)?/;(x)rfx = (ls j /s,o(x, u)ij (^X„+i(x, u) j dxrfu . 
Approximating R(X„) ^ R(X„+i), we obtain 



/s,o5V^/'(X„+i)(ix(iu 



Defining the implicit susceptibility x 



p(x)M(x)(I + R(x)) [I - N(x)] E„+i(x)rfx . 



^(^) = E + ^^'"W) [I - N(x)] p.(x) 



(47) 



(48) 



we have 

/ /5,o(x, u)5V(x, u)V^(X.„+i(x, u))(ix(iu = / V^(x)x(x)E„+i(x)(ix . 



(49) 



We treat the remaining integral by introducing the modified current 

j+(x)?/;(x)rfx = Qs J /s,o(x, u)V„+i/2 (x, u) tp (^X„+i(x, u) j rfxc/u . 

We then have 
q,At 



2eo 



8m^eo 
Sm^eo 



V X 

V X 



j+(x) X M(x) [I + R(x)] [I - N(x)] E„+i(x) V'(x)c/x 



7n+l/2(xj 



X [I + R(x)] [I - N(x)] 



E„+i(x) } ?/'(x)rfx 

(50) 



where use has been made of the identity U x U ® U = 0. We are then led to 
define the tensor ( as 



eo ^ rus 7n+i/2 



X [I + R(x)] [I - N(x)] 



(51) 



There follows 

2eo 



/,,o (5X ® V,+i/2 - V„+i/2 ® 5X) Vtjjdxdu = -At / V X (CE,+i)rfx. 

(52) 

Equation (i25D supplemented by Eqs. ( l46ll . (i49ll and (i52l) should be satisfied for 
any test function ip. As a result, we have to solve the local field equation 



En+l 



■ V X V X E„+i + xE„+i - AtV X (CE„+i) = Q , (53) 



where the source term now reads 

At- c^At^ 

Q = E„ j„+i/2 + c^AtV X B„„i/2 — V X V X E„_i . (54) 

eo I 

We have thus recovered the relativistic implicit method based on the Di scheme 
which was presented in Ref. [23|, with the only difference that the source term 
now involves the time- averaged field E„_i. It then appears that the direct implicit 
method^an be derived^ as a one-iteration Newton method with the starting values 
X(o) = X„+i, = U„+i/2 and E^") = 0. 
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3. Numerical resolution 



3.1. Resolution of the field equation 

In this section, we sketch the numerical procedure used to solve Eq. fl53ll in 
the case of a 2Dx-3Dv phase space with periodic boundary conditions along the 
transverse y axis. We have first to evaluate the implicit susceptibilities. These 
terms are computed for each macroparticle, yielding x(Xp,Up) and C(Xp,Up), 
before being projected onto the (x, y) grid through the usual formulas: 

^(^) = E E ^(^P - ' Up) > (55) 

s p 

^ W = E E ^(^P - ^)^(^^' Up) • (56) 



We then apply the iterative method of Concus and Golub [4j| to solve the elliptic 



system defined by Eq. (I53l) . which reads in the present case 

E^+i) + -^V X V X E('"+^) + xE(™+^) - AtV X (CE('"+^)) = Q^™) (57) 



The right-hand side of Eq. (1571 ) is given by 

qM = q _ _ + AtV X [(C - C°) E^"^)] (58) 

where m is the iteration index and x° and C° denote the y-averaged susceptibil- 
ities. The fast convergence of the scheme implies, in principle, slow variations 
of the field quantities in the y direction, but this has not proved particularly 
constraining for the physical situations we have considered. 

As is usual in electromagnetic PIC codes, two interleaved meshes are used for 
the spatial differencing of the grid quantities. The fields are discretized as fol- 
lows: Pij, Jz,i,j, Ez^iJ, Jx,i+l/2,j, Ex^i-\-l/2,j, By^i^i/2,j, Jy^ij+l/2, ,j+l/2 , -Bx,i,j+l/2 

and 5^,i+i/2,j+i/2- The x and C are stored at {i,j) except for Xii, Cii, C21, Csi, 
which are located at (z + 1/2, j), and X22, C12, C22, C32, located at + 1/2). 
Once space-discretized, the above equations are Fourier transformed along the y 
direction. Considering A^^^ grid cells, we obtain Ny one-dimensional equations to 
solve. Considering grid cells in the x direction, each equation gives a GN^ sys- 
tem of equations. These systems have a band-diagonal structure and are solved 
by a standard LU technique, using routines bandec and banbks of the numerical 



recipes library |45|. Details on spatial discretisations and Fourier transformations 



used to solve Eq. (I57l) are given in Appendix [Al 

3.2. Charge correction 

Our method to accumulate charge and current densities [Eqs. (|2T1) and (1471) ] 
does not satisfy charge conservation, which results into the violation of Poisson's 
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equation. This is a common flaw of early electromagnetic PIC codes [l| which 
may be corrected by a more sophisticated projection scheme |46|, l47|. A well- 
known alternative approach, which will be implemented here, is to correct the 
electrostatic part of the electric field E„+i solution of Eq. ( l53ll so that it fulfills 
Poisson's equation Using normalized quantities, our best statement of Gauss's 
law is 

V ■ E:+i = p„+i , (59) 

where E^+i represents the sought-for electric field. Using p„+i = Pn+i — V ■ 
(xE*^_^) , this can be reformulated as 

V ■ [(1 + x)K+i] = Pn+l ■ (60) 
Now, taking the divergence of Eq. fl53ll yields 

V ■ [(1 + x)E„+i] = V ■ Q (61) 

with generally V ■ Q 7^ Pn+i- We may first think of introducing a potential ip such 
that Q* = Q — fulfills V ■ Q* = pw+i, but this correction has been shown to 



cause spurious effects [20||. A proper correction makes use of the following form 

Q* = Q-(I + X)V^, (62) 

There follows 

V-[(l + x)V^] = V-Q-p.+i, (63) 

which is equivalent to 

V ■ [(1 + x) VV'] = V ■ [(1 + x)E„,+i] - pn+i , (64) 

where the only unknown is the scalar field ip. Eventually, the corrected field 
E*_^]^ ensuring Eq. ( l60l) is given by ^^+1 — En+i — Vip. Details on the numerical 
resolution of Eq. ( 1641) are given in Appendix [Bl 

3.3. Electromagnetic boundary conditions 

In this section we describe the implementation of injecting/outgoing boundary 
conditions on both sides of the simulation box. Incident and scattered electro- 
magnetic waves are assumed linearly polarized and depending on the phase term 
k • X — cut only. Waves polarized in the (a:, y) plane then verify 

Ef^ = 3^" cos Oi, (65) 

^scat ^ _QScat ^ ^gg^ 

where 6i and 6s denote respectively the incident and scattered angles. The total 
field becomes 



Ef = El^"' + Ef^ (67) 

rpinc 

cos 9i 



= * cos Os + (cos Oi + cos Os) (68) 
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Discretizing with centered finite differences in space and time gives 



4 {^y,l,j+l/2 + ^y,0,j+l/2 ^ ^y,l,j+l/2 + -^^1^,0^+1/2; " ^z,l/2,j+l/2 '^^^ 

+ ^-V2.>i/2 • (69) 

Using Maxwell-Faraday's equation, we can express -E'^q ^_,_^/2 ^ function of the 
field values at inner grid points and previous time steps. We have 

^s;>v2 = ^^;.t;+i/2 Os-i)-^-^ COS 9. (e-;2..i - ^:.Ia.) 

- 4^ COS ^.i?:7;£.,v2 + ^ (^s+i/2 - -+1/2 

- COS0, f f :7|, - f ."l + ^ (cos^, + cos^,) e:"<^'-+V2 



- ^ (Kl,i+l/2 + K\0J+l/2) , (70) 

where the coefficient A is given by 

A= fl + 2^cos^,) . (71) 



Ax 

A similar equation can be established for z-polarized waves, which reads 

2At 



^.,o,i-^^.,i,i^Aa;cose, 

45 „„„i/o 25At 



:iinc,n+l/2 / - COS ^.c 



where we have defined the coefficient B as 

f 2At , , 

V Ax cos 6',/ ^ ^ 

Note that the above equations only apply in vacuum. This is realized in 
practice by imposing boundary conditions on particles a few grid cells away from 
the outer boundaries of the computational domain. 



4. Numerical analysis of the adjustable-damping, direct implicit method 

4.1. Dispersion relation of electromagnetic waves in vacuum 

Our aim here is to quantify the error in phase velocity and the damping 
associated with electromagnetic waves as functions of the space and time steps. 
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In particular, we will demonstrate the possibility to control the wave damping by 
adjusting the parameter 6f. 

Combining Maxwell-Ampere's ( fT4l ) and Maxwell-Faraday's (fTSl) equations and 
assuming propagation in vacuum yield the wave equation 



En+i = 2E„ — E„_i — V X V X E„_,_i + E„_i 



(74) 



The time-filtered term involves the adjustable damping parameter 6f [Eq. ( ITTl) ] 
and can be expanded as 



E 



n+l 



E 



n-l 



E 



71+1 + " 



9f 



E 



n-l 



9f\ 6t 



E 



n-2 



(75) 



In a 2-D geometry, taking the electric field in the form E^ = Eo$(a;, t/)^;" with 
z = exp{—iuAt) and i = v^— T, Eq. (I75ll becomes 



En+i + E„_i = Eo$(x, y)<z ^ 



1 - 1 +^z + z^' 



OA 



1 - ^ 



l + ^z^^+ 1-^1 z-' + ... 



(76) 



where the adjustable damping parameter 9f G [0, 1]. Simplifying the series in the 
right-hand side of Eq. (1761) yields 



E„+i + E„_i = Eo<I>(x, y){z ^ 



iz'z + z^ 



1 - 



2 2.-^ 



/ 



(77) 



The electromagnetic wave is assumed polarized in the (x, y) plane with a harmonic 
dependence $(x, y) = exp ^{k^^ + ^yl/)]- Substituting Eq. (TfTl) into Eq. frf4ll and 
space-differencing the Laplacian leads, we get after some straightforward algrebra 
the following third degree polynomial equation 



= 2z 



1-— +^z + z' 



2 



e 



f 



2z-9 



f 



2 



^ , (78) 



where we have introduced 



Ax^ 



sm 



■ sm 



kyAy 



(79) 
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Equation (l78ll simplifies as 



\2 + -z^{A + Of) +z[2 + n\l- Of) + 29 f\ -9f = 0. (80) 
Let us first examine the special case 9f = 0. The roots of interest are solutions 

of 

z\2 + n^) -Az + {2 + n^) = (81) 

The discriminant A = 4 — (2 + fl'^Y being always negative, we get the roots 
z± = (2 ± iy/—A)/{2 + Q"^), which statisfy \z+\ = \z^\ = 1. We have there- 
fore demonstrated the absence of damping when 9f = 0. Figure [T] plots the 
normalized phase velocity = (where k = yZ/c^ + ky) for different values 
of cAt/Ax = cAt/Ay. The phase velocity error grows for increasing Ax and 
At/ Ax. A value cAt/Ax > 1, that is, violating the stability constraint of the 
standard explicit scheme, therefore implies a moderate spatial step k^Ax < 0.38 
(cAt/Ax = 1.27) so as to avoid excessive (> 5%) phase velocity error, which, in 



presence of relativistic particles, may cause unphysical Cerenkov radiation [48| . 

Let us now address the case of nonzero 9f. Figures [2] and [3] plot the normalized 
phase velocity v^/c (left) and damping rate Q^cuAt (right) of the least damped root 
of Eq. ( l80l) as functions of {k^Ax, kyAy) for ^/ = 1. Cuts of these two quantities 
in the plane ky = Q are represented in Figures [4] and [5] respectively. Again the 
phase velocity error grows for increasing Ax and At/ Ax. A value cAt/Ax > 1, 
therefore implies a reduced spatial step k^Ax < 0.28 (cAt/Ax = 1.27) so as to 
keep phase velocity error below 5%. In this case the damping rate, which also 
increases with Ax and At/ Ax, proves much too strong for applications relying 
on the propagation of an electromagnetic wave over several wavelengths. For 
example, assuming k^Ax = 0.2 and cAt/Ax = 1, a typical travel time of 200At 
requires IQ'a'IAt < 2.5 x 10~^ for a tolerable wave dissipation (< 5%). As seen 
in Fig. [5](right), this condition cannot be fulfilled when 9f = 1, which further 
demonstrates the need for an adjustable-damping scheme for a proper modeling 
of laser-plasma interaction. 

4.2. Dispersion relation of electrostatic plasma waves 

We will now focus on the numerical relation dispersion of the electron plasma 
fiuctuations in the case of a uniform, nonrelativistic Maxwellian plasma with a 
fixed neutralizing background. For this purpose, we shall adopt the formalism of 



Langdon [49i| that accounts for both finite space and time steps, as well as allows 
for an arbitrary time-differencing scheme of the Lorentz equation. An infinite 
number of macroparticles is assumed, yielding a continuous velocity distribution 
function (taken in the Maxwellian form). In this framework, as detailed in Ap- 
pendix O, the present adjustable-damping, direct implicit algorithm can be easily 
managed. The relation dispersion yielding the complex frequency a; as a function 
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Figure 1: Phase velocity of the least damped root of Eq. l(80l) as a function of {kxAx, kyAy), 
for different values of cAt/Ax cAt/Ay e {0.05, 0.66, 1.28, 1.9, 2.5} (from top to bottom) and 
9f = 0. A narrower {k^Ax, kyAy) range is represented on the right. 




Figure 2: Phase velocity (left) and damping rate QojAt (right) of the least damped root 
of Eq. l[80|) as a function of {k^- Ax , ky Ay) , for different values of cAt/Ax = cAt/Ay S 
{0.05,0.66,1.28,1.9,2.5} (from top to bottom on the left and bottom to top on the right) 
and Of = 1. 




Figure 3: Same as Fig. [2] but with a narrower (k^ Ax, kyAy) range. 
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1 2 3 1 2 3 

Ax kx A X 



Figure 4: Phase velocity (left) and damping rate 3a;Ai (right) of the least damped root 
of Eq. (|80|) as a function of (fcajAx), for different values of cAt/Ax = cAt/Ay G 
{0.05,0.66,1.28,1.9,2.5} (from top to bottom on the left and bottom to top on the right) 
and 9f — 1. Phase velocity without damping {df =0) is represented by dotted-dashed line. 




Figure 5: Same as Fig. |4]but with a narrower {kxAx) range. 
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of the wave number k then reads 



{kAx 



+ 



sm{kAx/2) 
kAx/2 

CJpAt)V2 



■E 

p=—oo 
+00 



sin {kpAx/2) 



2m+2 



sin(/cpAx) 



+00 



{kAx) 



siii(A:Az/2) 
kAx/2 



kpAx/2 

sin (/cp Ax/2) ^'"^^ sin(fcpAx 



5^ [1 + 



kpAx/2 



kpAx 



■S{9f) = 0, 



(82) 



where m is the order of the shape factor [l|. kp = k — 2ttp/Ax and ujq = 
uj — 2nq/ At are the ahased wave number and frequency, respectively. Z denotes 



the plasma dispersion function [50| whose argument is = Ug/ \f2kpVt (where 
is the electron thermal velocity). Moreover, we have defined the function S as 



m) = E 



s=0 



g-|(Ai3/Ax)2s2(fcAx)2(a;pAt)2 



(83) 



with the value 5(0) = 1. We have numerically solved Eq. (1821 ) using the nonlinear 



solver STRSCNE developed in Ref . [5ll| and the algorithm of Ref . [52| to compute 
the Z function. We will restrict the following analysis to systems characterized 
by a crude resolution of the Debye length (Ax/Xo > 1), as is commonplace in 
simulations of large-scale, high-density plasmas. 

Figure M displays the /c-dependence of the complex frequency of the fastest 
growing (or least damped) mode solution of Eq. ( !82i ) for 9f = 1, UpAt = 2 
and various values of Ax/Xd- For Ax/Xd = 32 {i.e., VfAt/Ax = 0.06), most 
of the /c-spectrum is damped except for a bounded unstable region located near 
kAx ~ 2.6 with a maximum growth rate ^ui/up ~ 0.011. This corresponds 
to the well-known finite-grid instability [l| commonly afflicting PIC simulations 
with Ax/ Ad ^ 1, and responsible for nonphysical field energy growth and plasma 
heating. This instability originates from the interplay of the aliased wave num- 
bers in Eq. (l82ll . Note also the nonphysical /c-dependence of the real frequency 
obtained at large uipAt : "Siuj is significantly below Up ai k = Q and further drops 
with increasing A;Ax. As seen in Fig. [6l decreasing Ax/Xd eventually leads to a 
complete stabilization of the system along with a displacement of the dominant 
mode towards low k values. For Ax/Xd = 4 {i.e., VtAt/Ax = 0.5), the least 
damped mode is thus located at kAx = 0.76 with ^uj/ujp ~ —0.1. This evolution 
points to a transition between spatial step-dominated and time-step-dominated 
regimes. 

The dependence of the characteristics of the dominant mode on the ratio 
Ax/Xo ^ 1 and the weight factor order is summarized in Table [H for Of = 1 and 
ujpAt = 2. The benefit of a high-order interpolation scheme is clearly evidenced: 
the system turns out to be entirely stabilized up to Ax/Xd = 32 with a quadratic 
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order 1 : MpAt - 2, &x/x^ - 32, 9 - 1 order 1 : WpAt - 2, Ax/lp - 20, e - 1 order 1 : WpAt - 2, Ax/J.p, -4,9-1 




Figure 6: Real frequency (blue) and growth rate (red) vs kAx of the dominant mode solving 
Eq. ^ with LUpAt = 2, Of 1 and a linear weight factor {n = 1): Ax/Xd = 32 (left), 20 
(center) and 4 (right). 



Ax/\d 


14.3 


22.6 


32 


64 


linear 


-0.024 


3.3 X 10-3 


0.011 


0.01 




(2.11) 


(2.42) 


(2.58) 


(2.85) 


quadratic 


-0.04 


-0.015 


-3.7 X 10-2 


2.8 X 10-3 




(1.96) 


(2.30) 


(2.48) 


(2.70) 


cubic 


-0.039 


-0.018 


-8.6 X 10-3 


-2 X 10-^ 




(1.84) 


(2.14) 


(2.36) 


(2.67) 



Table 1: Imaginary frequency '^uj/ujp (wavenumber kAx) of the dominant mode as a function 
of the ratio Ax/Xd and the weight factor order for ujpAt = 2 and 9f = 1. 



weight factor, and Ax/Xd = 64 with a cubic weight factor. In addition, the 
wavenumber of the increasingly damped dominant mode is shifted downward. 
A connection between the present calculations and previously published simu- 



lation results [l3|, l2l| is provided by Tables [2] and [3l which display the dependence 
of the dominant mode on the ratio VfAt/Ax = UpAt / {Ax / Xd) , as well as on the 
damping parameter (the time step being fixed to UpAt = 2). An extensive set 
of implicit electrostatic PIC simulations using the Di scheme {i.e., Of = 1) and 
linear interpolation has indeed revealed that satisfactory energy conservation can 
be achieved in the range [l3|, \21 



0-l<Vt^<l (84) 
Ax 

Even though the present stability analysis alone is not expected to account for 
the complex issue of numerical self-heating (l|,[5i|, the results of Table [2] are found 
in reasonable agreement with the lower bound of the above heuristic range, as 
they indicate a complete stabilization of the system for VfAt/Ax > 0.1 in case 
of a linear weigth factor and 9f = 1. For lower 9f values, stabilization is reached 
for increased VtAt/Ax. Moreover, Table [3] shows that the use of a quadratic 
weight factor permits to suppress the finite-grid instability at reduced VfAt/Ax 
(> 0.06 for 9f = 1). Similarly to Fig. [6l a clear transition from the high-/c spatial 
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n ni 1 
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(2.51) 


(2.51) 


(2.54) 


(^z.oo j 


0.1 


0.0204 


0.0185 


0.01 


-1.8 X 10-3 




(2.18) 

8 X 10"^ 


(2.18) 


(2.27) 


(2.33) 


0.25 


-7.4 X 10-3 


-0.04 


-0.08 




(1.05) 


(1.11) 


(1.28) 


(1.46) 


0.5 





-0.01 


-0.0508 


-0.105 




(0.39) 


(0.54) 


(0.63) 


(0.76) 


1 





-0.0102 


-0.0532 


-0.112 




(0.14) 


(0.27) 


(0.33) 


(0.39) 



Table 2: Imaginary frequency Sw/ojp (wave number /cAx) of the dominant mode as a function 
of the ratio wtAt/Ao; and the damping parameter Qj for WpAt — 2 and a hnear weight factor. 



regime to the low-/c temporal regime is evidenced when raising VtAt/Ax. As 
expected, a high-order (m > 1) weight factor, which enables to filter out high 
spatial frequencies, proves beneficial only in the high-/c, grid-instability regime 
(for VtAt/Ax < 0.25). Note that we have not considered values VtAt/Ax > 1 
since, in the present case, this would imply Ax/Xd < 2, a parameter range of 
little practical interest for the aforementioned applications. 

Further insight into the stability properties of the adjustable-damping scheme 
is given by fixing the ratio VtAt/Ax = 0.09 and varying accordingly the space 
and time steps. Equivalently, within the laser-plasma context which we propose 
to address, this can be achieved by fixing the parameters ujqAx/c and uoAt 
(where is the incident laser frequency) and varying the plasma density. The 
resulting data is displayed in Table [4] in the ranges 1.26 < Upt < 8.94 and 
14.3 < Ax/A/p < 101.1. One can see that a linear shape factor proves rather 
inappropriate for most of the parameter range considered. By contrast, complete 
stabilization is achieved for n > 2 weight factors. It is worth noting that, in 
terms of laser-plasma parameters, the rightmost column of Table [4] corresponds 
to a 2000nc, 1 keV plasma (where ric is the critical density at the laser frequency 
ujq) discretized with UoAt = 0.2 and uqAx/c = 0.1. In addition to accessing 
such extreme plasma conditions, employing a cubic weight factor may give the 
opportunity to reduce the damping parameter 9f. 
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0.1 
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VtAt/Ax 










0.05 


5.3 X 10"^ 


5 X 10-3 


3.5 X 10-3 


10-^ 




(2.54) 


(2.54) 


(2.58) 


(2.61) 


0.0625 


5.4 X 10-3 


4.8 X 10-3 


1.8 X 10-3 


-3.7 X 10-3 




(2.39) 


(2.39) 


(2.45) 


(2.48) 


0.1 


3.2 X 10^3 


1.1 X 10-3 


-8 X 10-3 


-0.0207 




(1.99) 


(2.02) 


(2.14) 


(2.24) 


0.25 





-8.1 X 10-3 


-0.039 


-0.078 




(0.81) 


(1.05) 


(1.22) 


(1.4) 


0.5 





-9.7 X 10-3 


-0.05 


-0.103 




(0.33) 


(0.54) 


(0.64) 


(0.76) 


1 





-0.01 


-0.053 


-0.11 




(0.14) 


(0.27) 


(0.33) 


(0.39) 



Table 3: Imaginary frequency Scj/wp (wave number klS.x) of the dominant mode as a function 
of the ratio vtAt/Ax and the damping parameter Of for ujpAt = 2 and a quadratic {n = 2) 
weight factor. 



Up At 


1.26 


2 


2.83 


3.46 


4 


5.66 


6.32 


8.94 


Ax/\d 


14.3 


22.6 


32 


39.1 


45.2 


64 


71.5 


101 


linear 


-0.0036 


0.0034 


0.0048 


0.0047 


0.0044 


0.0036 


0.0033 


0.0024 




(2.09) 


(2.41) 


(2.59) 


(2.67) 


(2.74) 


(2.85) 


(2.87) 


(2.96) 


quadratic 


-0.021 


-0.015 


-0.01 


-0.0078 


-0.0066 


-0.0044 


-0.0039 


-0.0026 




(1.95) 


(2.3) 


(2.5) 


(2.62) 


(2.68) 


(2.82) 


(2.85) 


(2.92) 


cubic 


-0.022 


-0.019 


-0.015 


-0.013 


-0.011 


-0.0079 


-0.0071 


-0.0051 




(1.83) 


(2.16) 


(2.36) 


(2.48) 


(2.56) 


(2.7) 


(2.76) 


(2.85) 



Table 4: Imaginary frequency Sw/wp (wave number kAx) of the dominant mode as a function 
of the space and time steps and the weight factor order, for a fixed ratio VfAt/Ax = 0.09 and 
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Figure 7: Propagation of a plane wave with 9f = 1 (top, left), 0f = (top, right), and a 
spatially varying 0f profile according to Eq. l[85|) (bottom). 



5. Numerical applications 

5.1. Wave propagation in vacuum 

Here, we illustrate the capability of the adjustable damping, implicit scheme 
implemented in the code ELIXIRS to manage the propagation of electromagnetic 
waves in vacuum. Let us consider a plane wave, with normalized vector potential 
ao = 3 and frequency uq, entering the left-hand side of a 1024 Ax x 4 Ay box, 
with Ax = 0.2c/uJo, Ay = O.Sc/cjq and At = 0.2c<Jq ^. The wave is injected and 
absorbed using the procedure detailed in 13.31 Figure [7](left) shows the expected 
monotonous damping of the incident wave induced when a spatially uniform 
damping parameter = 1 is applied. After propagating across the simulation 
box, the wave amplitude is measured to be 46% of the initial value, which is close 
to the theoretical value (49%). The opposite, dissipation-free case corresponding 
to ^/ = is displayed in Fig. [7](right). Finally, with the problem of laser plasma 
interaction in mind, we address the case of a spatially varying Of profile in the 
form 

ef = 0, 0<cuox/c<51.2 

9f = 1, 51.2 < iuox/c < 153.6 (85) 
9f = 0, 153.6 < uqx/c < 204.8 

Figure [7](center) shows that the discontinuity in does not cause significant 
spurious effects. This sought-for property is of major interest for modeling 
laser-plasma interaction as it allows the laser wave to travel unperturbed in 
vacuum over several wavelengths before reaching the overcritical target, whose 
numerical stability calls for finite numerical damping. For the sake of com- 
pleteness, we have checked that the weak (~ 0.1% in the present case) re- 
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flection arising at the discontinuity surface is consistent with Fresnel's formula 
R = (iV(i) _ iV(0))V(A^(l) + iV(0))2, where N{9f) = c/v^{9f) is the numerical 
refraction index derived in Sec. I4.1[ 

5.2. Plasma expansion into vacuum: benchmarking against explicit simulations 

As a first test of the implicit Vlasov-Maxwell solver, we simulate the dynamics 
of a plasma slab freely expanding into vacuum. The results of the implicit code 
ELIXIR are confronted to refined, explicit simulations performed with the code 
CALDER (sj]. We consider a O.Qc/ujp plasma slab composed of hot (10 keV) 
electrons and cold ions. In the implicit case, the simulation box is 103Ax x AAy 
large, with Ax = 2c/ Up and Ay = O.Ac/ujp (yielding the ratios Ax/Ad = 14 and 
VtAt/Ax = 0.14), whereas the explicit simulation handles a 1024Ax x 8Ay box, 
with Ax = Ay = 0.2c/ Up. A linear weight factor is used in all cases. 




Figure 8: Time evolution of the ion density profile: explicit (left) and implicit (right) simulations 
with Ax = 0.2c/ujp, At = O.lujp^, Np = 6 x 10^ and Ax = 2c/ujp, At = \ A^^ = 6 x 10'', 
respectively. The implicit damping parameter is 9f = 1. 
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Figure 10: Time evolution of the electron (red) and ion (green) kinetic energies: explicit (left) 
and implicit (right) simulations with Ax = 0.2c/a;p, At = O.lujp^, Np ^ 6 x 10^ and Aa; = 
2c/ujp, At = 2lJj7^, Np = 6 X 10'*, respectively. The implicit damping parameter is 6*/ = 1. 

Figures [U M and [10] plot the time evolution of the ion density profile, the ion 
phase space and the time evolution of the plasma kinetic energies, as simulated 
by the implicit and explicit codes. The implicit damping parameter is chosen to 
be = 1, whereas the total number of macroparticles Np is 6 x 10'' and 6 x 10^ 
in the implicit and explicit cases, respectively. Overall, albeit roughly resolved 
and strongly damped (as expected from Table [I]), the implicit scheme manages to 
satisfactorily capture the finely resolved, explicit results. Yet, the wave damping 
gives rise to artificial electron cooling, which results into a weakened ion accel- 
eration as seen in Figs. [9] and [TOl More quantitatively, the total energy drops 
by ~ 3%, yielding a maximum ion energy of ~ 160 keV, as compared to ~ 220 
keV in the explicit case. For the sake of completeness, we have carried out addi- 
tional calculations so as to assess the influence of the damping parameter and the 
number of macroparticules. For each simulation, we have measured the energy 
variation and the peak ion energy. The data thus obtained is summarized in 
Tables [5] and [6l The implicit scheme behaves reasonably well up to 9f = 0.15 
with an energy variation < 10%, comparable or better than its explicit counter- 
part for an equal number of macroparticles. Increasing the latter from 6 x 10^ to 
6 X 10^ approximately halves the energy variation but hardly changes the peak 
ion energy. The transition from numerical electron cooling and heating occurs 
between 9f = 1 and 9f = 0.5. Finally, the undamped (6'/=0) case is subject to 
a much stronger, if still limited, electron heating, which translates into a twofold 
overestimate of the peak ion energy. 

5.3. A parametric study of plasma self-heating and cooling 

We have carried out a series of simulations of the free evolution of an electron- 
ion plasma to gauge the potential discrepancy between the idealized linear anal- 
ysis of Sec. 14.21 and the actual predictor-corrector numerical scheme. Evidently, 
the objective is to gain further insight into the energy conservation properties of 
the latter and the predictive capability of the former. These calculations draw 
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AE/Eo Ion peak energy (keV) 



Explicit +9.3 % 232 

Implicit {Of = 1) -2.8% 162 

Implicit {9f = 0.5) +3.1% 208 

Implicit {9f = 0.15) +9% 273 

Implicit {9f = 0) +19.7% 451 



Table 5: Total energy variation and ion peak kinetic energy (keV) at 2600^ ^ with Np = 6x 10''. 



AE/Eq Ion peak energy (keV) 



Explicit +1 % 221 

Implicit {Of = 1) -1.4% 162 

Implicit {Of = 0.5) +1.5% 198 

Implicit {9f = 0.15) +4.5% 256 

Implicit {9f = 0) +12.4% 418 



Table 6: Total energy variation and ion peak kinetic energy (keV) at 2600^^ ^ with Np = 6x 10^. 



upon and extend the work of Ref. 2l|| to the electromagnetic regime. The 



system consists of a bounded electron-ion plasma with Tg = Tj = 1 keV and 
mi/irie = 900, extending over half a 300 Ax x A Ay simulation box. We have 
scanned the (Ax/A/), cjpAt) parameter space in the range [5,60] x [1,5]. In prac- 
tice, after introducing ljq, the frequency of a fictitious electromagnetic wave, 
and ric, the corresponding critical density, we have set Ax = 0.2c/ uq and var- 
ied the ratio rie/ric and the time step so that Ax/\d € {5,10,20,30,60} and 
(jjpAt e {1,2,5}. The damping parameter is Of = 1. The total simulation time 
is kept fixed at lOOOtUj^^. For each simulation, we have calculated the relative 
variation of the total kinetic energy per time step {AK / Kq) / N (where AK is the 
kinetic variation, Kq the initial kinetic energy and the number of time steps). 
To be complete, we have also performed electrostatic calculations, whereby the 
electric field is directly computed through the Poisson equation (l64ll . 
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60 


ojpAt 












1 


3.2 X 10-5 


3.2 X 10-^ 


1.1 X 10-3 


2.1 X 10-3 


4.7 X 10-3 


2 


-9.2 X 10-5 


1.5 X 10-^ 


7.9 X 10-^ 


1.5 X 10-3 


4.5 X 10-3 


5 





-1.6 X 10-^ 


1.2 X 10-^ 


4.7 X 10-^ 


1.7 X 10-3 



Table 7: Relative variation of the total kinetic energy per time step {AK / Kq) / N : electrostatic 
case and linear weight factor. 
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60 


OJpAt 












1 


2.8 X 10^5 


3.2 X 10-^ 


9.9 X 10-^ 


1.7 X 10-3 


2.8 X 10-3 


2 


-1.1 X 10"^ 


1.3 X 10-^ 


6.9 X 10"^ 


1.2 X 10-3 


2.6 X 10-3 


5 


-5.8 X 10-5 


-2.4 X 10"^ 


2.9 X 10^5 


3 X 10-^ 


9.5 X 10-^ 



Table 8: Relative variation of the total kinetic energy per time step {AK/Ko)/N: electromag- 
netic case and linear weight factor. 
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-3 X 10-5 


4 X 10-5 


2.3 X 10-^ 


4.6 X 10-^ 


1.1 X 10-3 
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-1.1 X 10-^ 


-3.5 X 10-5 


1.4 X 10-^ 


3.2 X 10-^ 


8.3 X 10-^ 
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-1.3 X 10-^ 


-2.2 X 10-^ 


-10-^ 





2.4 X 10-^ 



Table 9: Relative variation of the total kinetic energy per time step {AK/ Kq)/N: electromag- 
netic case and quadratic weight factor. 



The results are summarized in Tables (TH The associated plots of the kinetic 
energies are shown in Figs. [Tl]-[l3l each column corresponds to a specific value of 
Ax/Xd and each line to a specific value of UpAt. Note that we have excluded in 
these plots the case Ax/ Ad = 60 as it always gives rise to significant numerical 
heating. We have checked that the plasma kinetic energy makes up for most of the 
system energy. Overall, the electrostatic results prove close to the electromagnetic 
ones. Satisfactory energy conservation (< 10—*) is obtained for VfAt/Ax > 0.2 
and VfAt/Ax > 0.1 in the linear and quadratic interpolation cases, respectively. 
These lower bound values are in fairly good agreement, albeit slightly higher, 
with the linear results of Sec. 14. 2[ Larger ft At/ Ax ratios eventually lead to 
plasma cooling, 

5.4- High intensity laser interaction with an overdense plasma slab 
5.4- 1- Quasi- one- dimensional simulation 

Let us now address the problem of the interaction of a relativistic-intensity 
laser pulse with an overcritical plasma, which is the prime motivation behind this 
work. 

As a first illustration, we consider the case of a quasi-lD laser-plasma sys- 
tem. The irradiated target consists of a 60c/u;o-long, 1 keV, 200?t,c plasma slab 
preceded by a 18c/ci;o-long density ramp rising linearly from to 200nc . The 
incident electromagnetic plane wave has a 120ci;q ^ constant-intensity profile with 
a 22lJq^ rise time and a normalized amplitude = eEQ/mf^cojQ = 3. The 
implicit simulation employs a 2048Ax x AAy grid, with Ax = Ay = O.Ic/cuq 
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Figure 11: Time evolution of the total (blue), ion (red) and electron (green) energies: elec- 
trostatic case with linear weight factor. Ax/Xd = (5,10,20,30) from left to right and 
ujpAt ~ (1, 2, 5) from top to bottom. 
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Figure 12: Time evolution of the total (blue), ion (red) and electron (green) energies: elec- 
tromagnetic case with linear weight factor. Ax/Xd ~ (5,10,20,30) from left to right and 
ujpAt = (1,2, 5) from top to bottom. 
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Figure 13: Time evolution of the total (blue), ion (red) and electron (green) energies: electro- 
magnetic case with quadratic weight factor. 
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and At = O.McUq^, yielding, in terms of plasma parameters, Ax/Ad = 32 and 
ujp/S.t = 2 {vtAt/Ax = 0.06). The damping parameter in the electromagnetic 
solver, as well as in the particle pusher, is set to zero in the vacuum region and 
the moderately dense plasma region up to Ue = GOric, and to unity in the denser 
plasma region. Guided by the results of Sec. 15.31 we make use of a quadratic 
weight factor to reduce the numerical heating. The number of macroparticles 
per cell Np is varied from 100 to 1300. These calculations are compared with 
explicit simulations using the same parameters except for a decreased time step 
At = O.OSu;,^^ so as to fulfill the Courant stability condition. 







Explicit 


Implicit (Of = 0) 


Implicit {6 f = 1 if Tie > GOric) 


Np 


= 1300 


+14.4% 


+6% 


-3% 


Np 


= 400 


+15.3% 


+10.5% 


-1% 


Np 


= 100 


+22% 


+25.5% 


+12.7% 



Table 10: Quasi-ID laser-plasma interaction: energy variation in the explicit simulations with 
At = O.OBujq^ and the implicit simulations with At = O.lAuj^^ and varying 9f. See text for 
other simulation parameters. 



Table [TO] compares the values of the total energy variation (calculated after 
complete refiection of the laser pulse) as obtained in the explicit and implicit 
cases. Results from implicit simulations with zero damping are also displayed. 
Overall, except for Np = 100, for which case the three schemes behave similarly, 
the implicit simulations are found to achieve better energy conservation than 
their explicit counterparts. The benefit of a strongly damped scheme in the 
densest region of the plasma is mostly evidenced for Np = 1300 and 400. The 
not-so-good performances of the explicit calculations prompted us to carry out 
an additional, more refined explicit simulation that can serve more properly as 
a reference calculation. This simulation made use of a 4096Ax x 8Ay grid with 
Ax = Ay = 0.05c/u;o and At = 0.03a;Q ^, as well as of a third-order weight factor 
with Np = 650. It yielded a total energy variation of 4%. 

The electron {x,Px) phase space (integrated in the y-direction) is displayed 
in Fig. [14] for both explicit and implicit schemes. Consistently with the well- 
known ponderomotive heating mechanism arising at relativistic laser intensities, 
fast electrons are accelerated into the target as bunches separated by half the 



laser wavelength [55]|. The explicit simulation predicts maximum electron mo- 
menta about 20% higher than that predicted by the implicit simulation. Also, as 
a result of the damping of longitudinal beam-plasma modes, the implicit simula- 
tion exhibits a longer-lived separation between the thermal electrons and the fast 
electrons as the latter propagate through the target. In an actual solid-density 
configuration, though, the beam-plasma wave mixing observed in the explicit 
case should be suppressed by collisions as demonstrated in Ref. [H^. Yet, these 
discrepancies do not translate into major differences in the electron energy dis- 
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Figure 14: Electron {x,px) phase space at i = IQStJ,^^: explicit simulation (left) and implicit 
simulation with Of = 1 (right). In both cases, Np ~ 1300. See text for other simulation 
parameters. 




Figure 15: Electron energy distribution at different times: explicit simulation (red) and implicit 
simulation (blue). Energy is normalized by rrieC^. 
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Figure 16: Ion {x,px) phase space at t = 792luq ^: explicit simulation (left) and implicit simu- 
lation with Of = 1 (right). In both cases, Np = 1300. See text for other simulation parameters. 
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Figure 17: Time evolution of the electron (red) and ion (green) kinetic energies: explicit sim- 
ulation (left) and implicit simulation with 9f = 1 (right). In both cases, Np = 1300. See text 
for other simulation parameters. 

tribution as shown at three successive times in Fig. [151 In particular, the slope 
of the high-energy tail of the spectra is satisfactorily reproduced. The reduced 
electron heating gives rise in turn to a ~ 15% slower, space-charge-driven ion 
acceleration into vacuum as depicted by the ion {x,px) phase spaces of Fig. [161 

5.5. Two-dimensional simulations 

We now consider a fully two-dimensional laser-plasma system. The electron- 
ion plasma slab has a peak density of 200nc, a temperature of 1 keV and a 
thickness of 6c/ uq. A 12c/ci;o-long linear density ramp is added in front of the 
target. The simulation box consists of a 1024 x 512 grid with Ax = Ay = 
0.1c/ uq {Ax/\d = 32). The incoming laser pulse has unchanged parameters 
except for a 12c/u;o FWHM Gaussian transverse profile. Open and periodic 
boundary conditions are applied for the electromagnetic fields along the x- and 
y-axis, respectively. Due to memory constraints, we use a rather small number 
of macroparticles Np = 40. So as to stabilize the system, in addition to using 
a quadratic weight factor, the time step is significantly increased as compared 
to the previous simulations: At = 0.3u;o ^, which corresponds to cUpAt = 4.2 
and VtAt/Ax = 0.13. Particles are subject to periodic boundary conditions in 
the ^/-direction, and reinjected with their initial temperature in the x-direction. 
The damping parameter in the electromagnetic solver, as well as in the particle 
pusher, is set to zero in the vacuum region and the moderately dense plasma 
region up to = 30nc. Two maximum values of the spatially varying damping 
parameter have been tried in the denser plasma region: 6f = 0.1 and 0.5. The 
explicit simulation of reference makes use of a third-order weight factor with 
the parameters Aa; = Ay = O.OSc/cjq, At = 0.05ti;o^ and Np = 160. This 
parallel calculation takes 4.5h on 64 1.6 GHz Itanium 2 processors. By contrast, 
the (sequential) implicit simulations take 27h on a 2.66 GHz Intel Xeon X5355 
processor. 
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The time evolution of the particle kinetic energies is displayed in Fig. [181 All 
simulations predict about the same peak electron energy. Yet, the damped im- 
plicit calculations yield a faster decreasing electron energy. The total energy vari- 
ation, evaluated over the time interval 215 < uot < 715 (that is, after complete 
reflection of the laser pulse and before the fastest ions hit the box boundaries) 
is —12% and —15% for the 9f = 0.1 and 9f = 0.5 implicit cases, respectively, as 
compared to +5% in the explicit case. 




200 400 600 800 200 400 600 800 200 400 600 800 

ni„t 



Figure 18: Time evolution of the electron (red) and ion (green) kinetic energies: explicit simu- 
lation (left), implicit simulations with 9f ~ 0.1 (center) and 6*/ = 0.5 (right). 



Despite their crude time resolution and limited number of macroparticles, the 
implicit calculations manage to reproduce quite accurately the salient features 
of the fast electron and ion generation. This is evidenced by the electron and 
ion {x,px) phase spaces of Figs. fT9l and [20l as well as by the electron energy 
spectra of Fig. [22l As in the previous Section, if to a lesser extent due to the 
weaker numerical damping employed here, the implicit simulations somewhat un- 
derestimate the maximum electron energies. A 2-D picture of the fast electron 
generation is provided by the map of the electron kinetic energy density shown in 
Fig. [231 A reasonable agreement is observed between the three cases, each calcu- 
lation showing the characteristic 2u;o-bunched propagation of the fast electrons 
and their breakout into vacuum. 
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Figure 19: Electron {x,px) phase space at i = 96uJq^: explicit simulation (left) and implicit 
simulations with 9/ = 0.1 (center) and 9f = 0.5 (right). 
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Figure 20: Ion {x,px) phase space at t = 523^,^^: explicit simulation (left) and implicit simu- 
lations with 9f = 0.1 (center) and 9f = 0.5 (right). 
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Figure 21: Electron energy distribution at different times: expHcit simulation (red) and impHcit 
simulation with Of = 0.1 (blue). Energy is normahzed by nieC^. 
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Figure 22: Electron energy distribution at different times: expHcit simulation (red) and impHcit 
simulation with 9f — 0.5 (blue). Energy is normahzed by nieC^. 
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COqX/C {OqX/C 

Figure 23: Electron kinetic energy density (normalized by meC^ric) att = 67luq^ and t = 86ujq^: 
explicit simulation (top) and implicit simulations with 9f = 0.1 (center) and df = 0.5 (bottom). 
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6. Conclusion 



This paper has been devoted to the application of the relativistic direct im- 
plicit method to the problem of laser-plasma interaction. In contrast to closely 
related works 26|, 1271, l28|, our scheme, implemented inside the 2Dx-3Dv code 
ELIXIRS, allows for high-order weight functions and adjustable damping of the 
high-frequency waves. The latter capability, which extends to electromagnetic 
waves a method originally designed by Friedman [4ll| for electrostatic waves, per- 
mits to manage within a unified algorithm the dissipation-free, Courant condition- 
free propagation of the incident laser pulse through vacuum, while suppressing 
the need to resolve the high-frequency collective modes inside the dense plasma 
region. After having presented an original derivation of the adjustable-damping, 
direct implicit method as a simplified, one-iteration Newton scheme, we have 
carried out a thorough analysis of its numerical properties regarding both elec- 
tromagnetic and electrostatic waves. The latter study, accounting for the effects 
of finite At and Ax, the weight factor order and the damping parameter is found 
to provide useful hints when compared to the simulation results of the free evolu- 
tion of a plasma slab. Several numerical tests have been presented and successfuly 
benchmarked against finely resolved explicit simulations. In particular, we have 
demonstrated the ability of the code to capture the main features of the laser- 
plasma interaction despite cruder space-time resolution. Yet, our code being still 
sequential, its increased stability domain remains insufficient to access the large 
space- and time-scales managed nowadays by massively parallel explicit codes. 
The parallelization of our code is therefore required and will be the subject of a 
future work. 
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A. Numerical implementation of the field equation 

We detail here the numerical procedure to solve Eq. (157]) within a 2D ge- 
ometry. The Concus and Golub iterative method [4j| is applied to the three 



36 



components of Eq. (l57l) . The x-component writes 
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The ^-component writes 
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Assuming periodicity of the electric field along the y direction, we Fourier trans- 
form Eqs. (i86l)-(i88l) in this direction. We introduce E^ and El the real and 
imaginary parts of the Fourier transformed electric field. For notational simplic- 
ity, the index k will be omitted in the following. The real part of the Fourier 
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transform of Eq. (l86ll reads 
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The imaginary part of the Fourier transform of Eq. (l86l) reads 
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The real part of the Fourier transform of Eq. (l87l) reads 
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+ K")... {^C2f (-(^A,) + l)} + (Ei)^^^ {i^C2f sm(^A,)} = (Qf). 

(92) 
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The imaginary part of the Fourier transform of Eq. (l87l) reads 



+ (E^).-. {3^Ca°Bin(*A,)} + {-^Cf (co=(*A,) + l)} 

+ {--<kAy) - l) + 4^ (cos(^Ay) + l) - ^^1% (cos(A:Ay) + l) } 

+ (^i"). |-^sin(fcA|/)| + {Ei)^ 1^ (cos(fcAy) + l) | 

r p2A^2 _ 21,0 ^ A . 

+ (^3.^1/2 |^A^A^^i^(^^^) - V^^"^^^^) - ^Cl%sin(fcA^)| 

+ - + 4^ + + ^Cv2 (cos(fcA,) + l) } 

l^Wi+i \ 2Aa;2 2Aa:^*+i j 
+ K")..! sm(^A,)} + (Ef),^^^ {^C2f (cos(^A,) + l)} = (Qi). 

(93) 
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The real part of the Fourier transform of Eq. (l88l) reads 
, / , 23,0 1 

r ~ Ay- ~ 1 
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The imaginary part of the Fourier transform of Eq. (l88l) reads 



/ T-i/X J Cp' 23 



+ (S,'). 1^ (co8(*Ai/) + l) + ^C"'° (cos(fcAj,) - l) I 



+ (^^f).{^e°=ta(*A!,)' 

+ (^.')«{-^-^cS°}-(«0, (») 

Considering A'^. grid points along x-direction Eqs. (I90l)-(l95l) can be formulated as 
a band-diagonal system of equations, which we solve using a LU technique [45 
for each of the Ny modes of the discrete Fourier transform. Then we compute 
the field solution in real space by inverse Fourier transformation. 

B. Numerical implementation of the charge correction step 



We detail here the numerical procedure to solve Eq. (1641 ) within a 2D geom- 
etry. As for the wave equation, we make use of the Concus and Golub iterative 



method [4j|, which writes in the present case 



- V ■ [(1 + /)VV^("^+^)] = p - V ■ [(1 + x)E„+i] + V ■ [(x - X°)V^(")] (96) 

where x° = [x'^''^] ;<3 denotes the y-averaged x susceptibility tensor with 
^ki,o _^ ^ki E„+i is the solution of the wave equation (l53l) and m denotes 
the iteration index. Omitting the latter, we discretize the above equation in the 
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form 



Ax 
1 

2 Ax 
1 

'2K~y 
1 



1 + X 



11,0 



12,0 1 / ; 



21,0 



'2 Ay 
1 

X 



11,0 



^■+'2A 



i+l,i+l 



12,0 

Xi-i,j 



; N 21,0 

- ^i-i,j+i) - Xi,j-i 



Xi-l/2,jJ 
1 



2Ay 
1 

2Ax 



22,0 
^j,i-l/2 



(V^i-lj+l - V'i-l,j-l 



(97) 



where we have defined the source term 



+dy [ix'' - x'''')d.i^ + ix'' - x'''')dy^] + P 

-d, [(1 + x'')E.] - ix'^Ey) - {x''E,) 



(98) 
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A centered spatial discretization of Eq. ([98|1 is given by 



Si,j — + * 



+ 



1 

Ax 
1 

2Ax 



Xi 

1 



12,0\ 



12 _ 

i+i,j ^m/2Ay 



) 



2Ay 



+ 



1 

1 

Ax 
1 

1 

2Ax 
1 



[(1 + Xi+l/2,j) -^1^,^+1/2 J - (1 + Xi-l/2,i)-^2^,i-l/2,j] 



A«+lj 



,12 



(-£^S/,j+l,i+l/2 + Ey^i+l,j~-l/2) (-Ej/,i-l,j+l/2 + -Ey,i-lj-l/2) 



2Ay 
1 

A^ 
1 



Aij + l 



x,i+l/2,j+l + -E'x,fi-1/2J+1 







Ajj-l 



2A?/ 



[(l + xS+1/2) ^y,i,j + l/2 - (1 + x3-l/2) Ey^iJ-1/2] 

[Xij+iEz,ij+i - x'ij^iEz,i,j-i] 



{E, ,j+l/2J-l + -E'a;,i-l/2,j-l) 



(99) 



The above equations are Fourier transformed along the y direction. Considering 
Ny grid cells we have to solve Ny one-dimensional equations. Assuming grid 
cells in the x direction, each equation turns out into a 2Nrc system of equations. 
These systems have a band-diagonal structure and are solved with a LU technique 



45 , 



C. Derivation of the dispersion relation of electron plasma waves w^ith 
finite Ax and At 

We restrict our analysis to a one-dimensional, nonrelativistic electrostatic 
plasma with immobile ions. In the following, we adopt the methodology and 
notations of Ref. For a single macro-particle, the adjustable-damping scheme 
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(l7ll- (fT0l) can be formulated as 



22 



23 



an-3 + 



(100) 



where n stands for the time step index. We now assume a harmonic form for the 
interpolated electric force F^^^ = F(fc)e*^^^~'^*). As a direct consequence of the 
PIC interpolation scheme, we have the relation [H 

F{k) = qE{k)S{-k) (101) 

where E{k) and S{k) are the discrete Fourier transforms of the electric field and 
the m-order weight function, respectively. The latter reads 



S{k) 



" sin(fcAx/2) l '"^^ 
kAx/2 

The first-order acceleration term can then be expressed as 

m 



m 

m 

m 

m 

m 



exp [i{kx^J^^ — Lutr 



exp [ik{xo + v^^hn) — iootn] 
exp ikxo exp [i{kv — uj)nAt] 



Defining A{k) = ^(^^^^ 
_ At^ 

2Xj^ -|- Xn—1 



gifcxo ^ ^ ^i(kv-uj)At^ ([T0311 reads 

2 



2 
At2 



A{k) < 2"+^ + -z" + 



A{k)z'''' z-^ 



'f 2 

-^z + Z^ 



+ 1- 



/ \ .-2 



1 + ^^-1 + 



This equation can be further simplified as 

At^ 



A{k)2z 



2z-e 



(102) 



(103) 



(104) 



(105) 
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We linearize Xn = Xn^ + Xn^ where Xn^ = Xg^^ + vll^hn 



Where the polynomial V reads 



A{k)V{k) 



V{k) = 2z' 



2z-ef 



(106) 

(107) 
(108) 



We deduce that Xn\xo.,vo,tn) varies as e**^^^ uj)nAt _ Hence we find the 
solution 



X. 



(1) ^ ^±!_p(J^yi{kx-ujt) 



m 



(109) 



To evaluate the charge density, we introduce the dipole density 
P{x,t) = noq / dvfo{v)x^^\x,v,t) 



noq 



m 



F{k)e 



i{kx—LL>t) 



dvf oiv 



2m 



(Asin(^-Hf)' 



(110) 



The first and second terms of the right-hand side correspond to the explicit 
leapfrog scheme and the implicit correction, respectively. Assuming a Maxwellian 

1 



distribution /o(f) 



Vt\/2lT 



exp 



ft 



the latter can be written 



s=0 



(2/0,)' ^ die,) 



E 



E 



t my 



dvfo{v)e-'''^'^' 
J^{fo){ksAt) 



E 



-e 2 



f (sA;At)2 



where denotes the Fourier transform. Thus the polarisation becomes 



(111) 



P(^,t) = !^F(fc)e*(*^^--*) — 
m 4 



/^(.)^cotan 



(uj-kv)^ 



dv 



°° AujsAt 2 



2m 



^ (2/^/) 



-e 2 



(112) 
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We can develop cotan as a series in the form 
cot an 



{uj-kv} — 



— y 



At ^ uj — kv — qijja 



(113) 



The continuous charge density is given by pp = —V ■ P, which writes in Fourier 
space Pp(A;) = —ikPik). The discrete charge density is then given by 

V 

— ^ ^ ^ kpS {kp) P (^kp) 
p 

-^y\sikp)\^^Eikp) y [dv- 



— t 



u — kpV — qUg 



°o iujsAt „2 

e -^(sfcAt)2 



s=0 



{2/6) 



-e 2 



(114) 



Using centered space-differencing, discrete Fourier transform of the relation E 
—d(j)/dx gives 

E{k) = -iK{k)(t){k) = -iK{k)(t){k) , 

where 

sin(/i;Ax) 



K{k) = k- 



kAx 



(115) 
(116) 



The Poisson equation as modified by the direct implicit method reads 



V ■ (V0„H 



Pn+l 



Centered space-differencing followed by a Fourier transformation gives 



K\k)^{k) 



eo 



where we have defined 



K\k) = e 



sin(A;Ax/2) 



1 2 



(117) 



(118) 



(119) 



kAx/2 

Combining Eqs. (1114p -( fTT9l) . we obtain the dispersion relation for an infinite elec- 
trostatic one dimensional plasma taking into account both spatial and temporal 
discretizations 

,2 



€{uj,k) = 1 + 



+ 



UJr, 



K^{k) 



y\S{kp)\'K{kp) J2 J 



dfo{v)/dv 



q=~OD 



kpV — qujg 



ui Ae 

K'^ik) 2 



ykp\S{kp)\'Kikp)yj^e--^^'^^''^'^' = 0, (120) 



s=0 
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where kg = 211/ Ax, ujg = 27r/At, Ug = u — qug and hp = k — pkg. 
Exploiting the Maxwellian form of /o, we have 

dv^^^ = [1 + UZ{Q] , (121) 

UJq — kpV kpVf 

where C,q = ^ and Z denotes the Fried and Conte plasma dispersion function 
Z [50|, defined by 

du^- with $^(^) > . (122) 
00 ~ s 

Finally, substituting Eq. ffT2TD into Eq. [I20] yields Eq. fl82|) . 
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